// Trabalho Final de Programacao Concorrente
// Jacobi-Richardson
// Daniel Hatanaka, 7277579
// Guilherme Kalatzis, 7239426
// Thiago Tanaka, 7153014
// Prof. Dr. Julio Estrella, Estagiario PAE Luiz Henrique Nunes

//Entrada de dados segue
//3 -> ordem da matriz (J_ORDER)
//2 -> fila para ser avaliada (J_ROW_TEST)
//0.001 -> erro permitido (J_ERROR)
//20000 -> número máximo de iterações (J_ITE_MAX)
//4 2 1 -|
//1 3 1 -|-> matriz A (MA)
//2 3 6 -|
//7   -|
//-8  -|-> matriz B (MB)
//6   -|
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

//Definicao do numero de threads
#define NTHREADS 8

int main(int argc,char *argv[]){
		//Setando o numero inicial de threads
		omp_set_num_threads(NTHREADS);
		//Declaracao de variaveis, matrizes e tempos.
        float *AUX;
        float *AUX1;
        float **MA;
        float *MB;
        float tempo;
        double inicial = omp_get_wtime();
        double tempoAloc;
        double inicialAloc;
        double tempoInsercao;
        double inicialInsercao;
		float erro;
        float dif;
        float numerador = 0.00 ;
        float denominador = 0.00;
        float max = 0.00;
        float J_ERROR;
        float result;
        int J_ORDER;
		int J_ITE_MAX; 
		int J_ROW_TEST;
        int i;
        int j;
        int count = 0;
        int flag = 0;
		//Declaracao do arquivo de entrada
        FILE *pFile;
		char nome_arquivo[15];
        printf("Nome do arquivo: ");
        scanf("%s", nome_arquivo);

        pFile = fopen ( nome_arquivo , "rb" );
        if (pFile==NULL) {
        	fputs ("File error",stderr); exit (1);
        }
        //Voltamos o arquivo para o inicio
        rewind (pFile);
        //Leitura das variaveis
        fscanf(pFile, "%d", &J_ORDER);
        fscanf(pFile, "%d", &J_ROW_TEST);
        fscanf(pFile, "%f", &J_ERROR);
        fscanf(pFile, "%d", &J_ITE_MAX);
		inicialAloc = omp_get_wtime();
		//Alocacao de matrizes e vetor B
		//Malloc das variaveis
        MA = (float **)malloc(sizeof(float *)*J_ORDER);
        for(i = 0;i < J_ORDER;i++){
                MA[i] = (float *)malloc(sizeof(float)*J_ORDER);
        }
        MB = (float *)malloc(sizeof(float)*J_ORDER); 
        //Tempo final é o tempo inicial de alocacao menos o final                          
        tempoAloc = omp_get_wtime() - inicialAloc;
        printf("----------------------------------------------------------------\n");
        printf("Tempo de alocacao: %.10lf seg\n", tempoAloc);
        printf("----------------------------------------------------------------\n");
        //Tempo inicial de insercao de dados
        inicialInsercao = omp_get_wtime();

        for(i = 0;i < J_ORDER;i++){ 
                for(j = 0;j < J_ORDER;j++)
                		fscanf(pFile, "%f", &MA[i][j]);
        }

        for(i = 0;i < J_ORDER;i++){
                fscanf(pFile, "%f", &MB[i]);
		}
		//Tempo final
        tempoInsercao = omp_get_wtime() - inicialAloc;
        printf("Tempo de insercao de dados: %.10lf seg\n", tempoInsercao);
        printf("----------------------------------------------------------------\n");
        inicial = omp_get_wtime();
		//Malloc dos vetores AUX e AUX1
		for(i = 0;i < J_ORDER;i++){
                AUX = (float *)malloc(sizeof(float)*J_ORDER);
		}
        for(i = 0;i < J_ORDER;i++){
                AUX1 = (float *)malloc(sizeof(float)*J_ORDER);
				AUX[i] = 0;
		}   
		i=0;
		j=0; 
 
        
        
		//O critério de parada:
		//Atingir o erro (J_ERROR) ou
		//Atingir o número de iterações máximo (J_ITE_MAX)
		
		//Metodo Jacobi Richardson
        do{
        	//Denominador e numerador do calculo do erro
            denominador = 0.00;
            numerador = 0.00;
            count = count + 1;
            //loop principal do jacobi-richardson
            for(i = 0;i < J_ORDER;i++){
                    AUX1[i] = 0.0;
					
                    for(j = 0;j < J_ORDER;j++){
                             if(i != j){
                                 AUX1[i] = AUX1[i] + MA[i][j] * AUX[j];
                            }
                     }
                    AUX1[i] = 1 / MA[i][i] * (MB[i] - AUX1[i]);
                    dif = AUX1[i] - AUX[i];
                    if(dif < 0){
                        dif = dif * (-1);
                    }
                    if(numerador < dif){
                        numerador = dif;
                    }

                	dif = AUX1[i];

                    if(dif < 0){
                        dif = dif * (-1);
                    }

                    if(denominador < dif){
                            denominador = dif;
                    }
            }
            for(i = 0;i < J_ORDER;i++){
                    AUX[i] = AUX1[i];
			}
            erro = numerador/denominador;

        }while(erro >= J_ERROR && count < J_ITE_MAX);
				//Paralelizacao do calculo do resultado
				#pragma omp parallel for
                for(i = 0; i < J_ORDER; i++){
                	#pragma omp critical
             		result =  result + MA[J_ROW_TEST][i]*AUX1[i];

        		}
        printf("----------------------------------------------------------------\n");
        printf("Iterations: %d\n", count);
        printf("RowTest: %d => [%f] =? %f\n", J_ROW_TEST, result, MB[J_ROW_TEST]);
        printf("----------------------------------------------------------------\n");
		tempo = omp_get_wtime() - inicial;
        printf("Tempo de execução do Jacobi-Richardson: %.10lf seg\n", tempo);
        printf("----------------------------------------------------------------\n");
        //Libera alocacao de dados
        free(MA);
        free(MB);
        free(AUX);
        free(AUX1);
        printf("---------------------------Fim Programa-------------------------\n");
        return 0;
}

